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ABSTRACT 

The spectral slope of strong MHD turbulence has recently been a matter of controversy. While 
Goldrcich-Sridhar model (1995) predicts Kolmogorov's -5/3 slope of turbulence, shallower slopes were 
often reported by numerical studies. We argue that earlier numerics was affected by driving due to 
a diffuse locality of energy transfer in MHD case. Our highest-resolution simulation (3072 2 x 1024) 
has been able to reach the asymptotic —5/3 regime of the energy slope. Additionally, we found that 
so-called dynamic alignment, proposed in the model with -3/2 slope, saturates and therefore can not 
affect asymptotic slope. The observation of the asymptotic regime allowed us to measure Kolmogorov 
constant Cka = 3.2 ± 0.2 for purely Alfvenic turbulence and Ck = 4.1 ± 0.3 for full MHD turbulence. 
These values are much higher than the hydrodynamic value of 1.64. The larger value of Kolmogorov 
constant is an indication of a fairly inefficient energy transfer and, as we show in this Letter, is in 
theoretical agreement with our observation of diffuse locality. We also explain what has been missing 
in numerical studies that reported shallower slopes. 



1. INTRODUCTION 

The equations of incompressible ideal magnetohydro- 
dynamics, written in terms of Elsasser variables, 



9 4 w ± + 5(w^ 



V)w ± = 0, 



(1) 



where S = (1 — VA 1 V) is a solenoidal projec- 
tion, and w ± (Elsasser variables) are w + = v + b and 
w~ = v b, b = B/(47rp) 1 / 2 are remarkably similar to 
the Euler's equation. However, have different trans- 
formation properties than v. While in hydrodynamic 
turbulence the local average velocity can always be ex- 
cluded by the choice of reference frame, in MHD al- 
ways contain the average local mean magnetic field that 
can not be excluded. This leads to a situation when 
large scale magnetic field is much stronger than small- 
scale turbulent perturbations and dynamic s is domi- 
nated by this local m ean magnetic field (jlroshnikovl 
119631 : iKraichnanl [1965) . But turbulence does not be- 
come weaker down the cascade as was proposed in 
afor ementioned models . A proper perturbation the- 
ory (jGaltier et all 120001) revealed that MHD turbulence 
has a tendency of becoming stronger on smaller scales, 
rather than weaker, due to the fact that the cascade in- 
creases perpendicular wavenumber kj_, keeping parallel 
wavenumber ku constant. The "strength" of turbulence, 
as the ratio of the mean-field term to the nonlinear term 
can be approximated as £ = wkj_/vAku and can increase 
due to the increasing anisotropy of perturbations k±/kn 
down the cascade. As turbulence becomes marginally 
strong (£ ~ 1), i.e., the linear term is comparable to the 
nonlinear term, the cascading timescales become close to 
the dynamical timescal es r ease ~ Tdvn = l/wk±. How- 
ever, as was arg ued in iGoldreich fc Sridharl (|1995D . the 
perturbation frequency lo has a lower bound due to an 
uncertainty relation □ r casc Lu > 1, therefore the combi- 

1 Another bound on u> t hat follow from the uncertai nty in the 
direction of the vector ( Beresnyak & Lazarian 2008 ) , give the 
same estimate in the case of balanced turbulence considered in 



nation of this lower bound, that limits the strength of 
turbulence £ and the tendency of turbulence to become 
stronger will make it "critically balanced" with £ ~ 1. 
This critically balanced cascade is strong in a sense that 
cascading timescale is always of the order of the dynamic 
timescale and will, therefore, have a Kolmogorov -5/3 
spectrum. These critically balanced perturbations will 
be strongly anisotropic with respect to the local mean 

— 1/3 

magnetic field. Using £ ~ 1 and w ~ k ± ' one obtains 

fen ~ fc 2 ^ 3 , i.e. the anisotropy will increase towards small 
scales without limit. One can further simplify Eq. 1 by 
neglecting the term (Sw^V^Sw^ which is much smaller 

than the mean field term (v^V||)(5k; ± . After this Eq. Q] 

splits into two equations, one for 5w^, which, in this 

strongly anisotropic case, fey <C k± represents slow (or 

pseudo-Alfven) mode and the equation for Sw^ which 
represent Alfvenic mode. The equation for slow mode is 
passive and does not provide any back-reaction for the 
Alfvenic equation which could be written in the following 
form: 



d t 5wf T (v A • V|| )Sw^ + S(Swf • V_l)5w± = 0, (2) 

One can study this purely Alfvenic dynamics and 
assume that the omitted slow mode has similar cas- 
cade and similar statistical properties. Eq. [2] is known 
as reduced MHD app roximation or R MHD, see, e.g., 
iKadomtsev fc Pogutsd (fl97l : iStraussI (fl976l ). RMHD 
equations provide further support towards Goldreich- 
Sridhar model, as it has a precise symmetry with respect 
to anisotropy and the strength of the mean field. Indeed, 
as long as one increases va and stretches the fields in 
the parallel direction, decreasing V||, by the same factor, 
Eq. 2 will be unchanged. Furthermore, a Kolmogorov 
argument of universality of nonlinear dynamics at each 



Goldrcich & Sridhar (1995). In the imbalanced case it could lead 
to a modified relation. 
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TABLE 1 
Three-dimensional simulations 



Run 


Type 


T^x ' T^y ' T^z 


Dissipation 


<^> 


L/ V 


HI 


hydro 


512 3 


-3.02 ■ 10~ 4 fc 2 


0.091 


190 


H2 


hydro 


1024 3 


-1.20 ■ 10~ 4 fc 2 


0.091 


370 


Ml 


MHD 


1024 3 


-1.63 • 10" 9 fc 4 


0.159 


280 


Rl 


RMHD 


256 ■ 768 2 


-6.82 ■ 10~ 14 £i 6 


0.073 


280 


R2 


RMHD 


512 ■ 1536 2 


-1.51 ■ 10~ 15 fc 6 


0.073 


570 


R2.5 


RMHD 


1536 3 


-1.51 ■ 10- 15 fc 6 


0.073 


570 


R3 


RMHD 


1024 ■ 3072 2 


-3.33 ■ 10" 17 fc 6 


0.073 


1100 



scale, which is based on a two-parametric scaling sym- 
metry, could be amended it with a proper scaling for the 
anisotropy: 

w^wi, A — > XB, t^tB/A, A — s> AB/A, (3) 

where A is a perpendicular scale, A is a parallel scale, 
A and B are arbitrary parameters. Due to this precise 
symmetry one can hypothesize that strong Alfvenic tur- 
bulence has a universal regime, utilizing the same argu- 
mentation as IKolmogorovl (|T941). In nature, this univer- 
sal regime can only be achieved as long as 5w ± <C va- 
In numerical simulations, we can directly solve the re- 
duced Eq. [2] , which has precise symmetry already built 
in. From practical viewpoint, the statistics from the full 
MHD simulation with 5w ~ O.Iva is virtually indis- 
tinguishable from RMHD statisti cs and even Sw^ ~ va 
are fa irly similar to the former (jBeresnvak fc Lazarianl 
I2009aj) . Note, that as we see above, both MHD and 
RMHD dynamics are essentially three-dimensional. In 
this paper we use both full MHD simulations and RMHD 
simulations. Statistically isotropic MHD simulation is 
used to determine a fraction of total energy contained 
in the slow mode, while RMHD simulations are used to 
study properties of the universal Alfvenic turbulence. 

Previous numerical work confirmed scale-dependent 
anisotropy of the strong MHD turbulence (se e , e.g., 
ICho fc Vishmacl l200fl iMaron fc Goldreichl l200l . The 
precise value of the spectral sl ope, however, was a mat- 
ter of debate. In particular, (jMuller fc Grappinl 120051 ) 
claimed that the mean field strong turbulence has a 
slope of —3/2. This m otivated adjustments to the 
Gold r eich-Sridhar mode l (|Galtier et all 120051: iBoldvrevi 
120051: iGogoberi dzc 2007 3). A model wit h so called "dy- 
namic alignment" (Boldyrev 20Q5|, 120061 ) became popular 
after the alignment was discovered in numerical simula- 
tions (Beresnv ak fc Laza rian 2006). Boldyrev model as- 
sumes that the alignment between velocity and magnetic 
perturbations decreases the strength of the interaction, 
also it assumes that the alignment is a power-law func- 
tion of scale, increasing indefinitely towards small scales, 
modifying the spectral slope of MHD turbulence from 
the —5/3 Kolmogorov slope to —3/2 slope. In this pa- 
per we debate the assumption that the alignment is a 
power-law function of scale. We also prove that earlier 
measurements of the slope were premature and were un- 
able to reach asymptotic slope due to diffuse locality of 
MHD turbulence. 

2. NUMERICAL METHODS 

We used a pseudospectral dealiased code that was 
able to solve hydrodynamic, MHD and RMHD equa- 
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Fig. 1. — Our hydrodynamic simulations reproduce C'k = 1.64. 
As long as the turbulence is local and the effects of large scales 
could be neglected in the inertial range, the larger 1024 cube could 
be seen as consisting of a eight smaller 512 cubes, as on the right. 

tions. The RHS of Eq. 2 has an explicit dissipation 
term — u„(— V 2 )" / ' 2 w ± and forcing term f. The code 
and the choice for numerical resolution, driving, etc, 
was described in great detail i n our ea rlier publications 
(jBeresnvak fc Lazarianl l2009allE f2010h . Table 1 shows 
the parameters of the simulations. The Kolmogorov 
scale is defined as 77 — (V^/e) 1 ^ 3 ™" 2 ), the integral scale 
L = 3ir/4E J Q k~ 1 E(k) dk (which was approximately 
0.79 for Rl-3). Dimensionless ratio L/rj could serve as a 
"length of the spectrum" , although usually spectrum is 
around an order of magnitude shorter. 

The resolution in the direction parallel to the mean 
magnetic field, n x was reduced by a factor of 3 for sim- 
ulations Rl-3. This was possible due to an empirically 
known lack of energy in the parallel direction in fc-space. 
We ran a simulation R2.5 which has full resolution in 
n x to compare with R2 and check the influence of this 
resolution reduction on the power spectrum. Although 
the bottleneck effect was slightly less pronounced in R2.5 
compared to R2, there was only a small influence in the 
inertial range. We concluded that using n x reduced by a 
factor of 2 or 3 is possible. 

For the purpose of this paper we used driving that had 
a constant energy injection rate. In RMHD simulations 
Rl-3 we drove turbulence to the amplitude that it will be 
strong on the outer scale. Rl-3 were started from lower- 
resolution simulation that reached stationary state and 
were further evolved in high resolution for approximately 
12 Alfvenic times, which, for strong MHD turbulence also 
correspond to about 12 dynamical times. The averaged 
quantities were obtained for the last 6 Alfvenic times. 
In all magnetic simulations Ml, Rl-3, we were using hy- 
perviscosity (n > 2) instead of normal viscosity. This is 
possible due to the fact that bottleneck effect is much 
less pronounced in the MHD case, compared to hydro. 

3. SPECTRA AND UNIVERSALITY 

Much of the study of hydrodynamic turbulence was 
dedicated to Kolmogorov model whic h assumes a univer - 
sal cascade of energy through scales (Kolmog orovlll94l|) . 
This model predicts that the power spectrum of turbu- 
lence, E(k), will be a power-law function of scale, 

E(k) = C K e 2 / 3 k' 5 / 3 . (4) 

where Ck is a Kolmogorov constant. It is well-known 
that this scaling is not precisely correct and typically 
has an intermittency correction (kL) a , where L is an 
outer scale and a i s a small number, around 0.035 
IjShe fc Leveaueiri994f ). However, in simulations or mea- 
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Fig. 2. — Resolution study for MHD simulations. As MHD is 
less local than hydro, convergence require higher resolution. The 
estimate of Kolmogorov constant from averaged spectra for this 
purely Alfvenic turbulence is Cka = 3.2 ± 0.2. 

surements with small inertial range this correction can 
often be neglected. In particular, a compilati on of experi- 
ment al results for hydrodynamic turbulence ()Sreenivasanl 
1995) suggests that a Kolmogorov constant is univer- 
sal for a wide variety of flows. High-resolution numeri- 
cal simulati ons of isotropic inc ompressible hydrodynamic 
turbulence (|Gotoh et alll2002l ) suggest the same value for 
the Kolmogorov constant. 

A robust method for determining the spectral slope 
and the Kolmogorov co nstant from simul ations is a res- 
olution study (see, e.g.. lGotoh et all l2002'). when a num- 
ber of numerical experiments are performed with differ- 
ent resolution and the spectra are plotted with respect to 
the dimensionless wavevector, kr\. A physical meaning of 
such a comparison is based on an assumption that a simu- 
lation with higher numerical resolution can be considered 
both as a simulation resolving smaller physical scales and 
as a simulation of a larger volume of turbulence (see Fig. 
[T|). This assumption is true as long as turbulence can 
be considered local, i.e. the effects of driving can be 
neglected in the inertial range. Our hydrodynamic simu- 
lations reveal a good convergence of spectra with numeri- 
cal resolution and show a universal Kolmogorov constan t 
consistent with the one obtained in iGotoh et all ([2002). 
Also the shape of the dissipation range is similar to the 
one in aforementioned paper, showing a typical "bump" 
due to a bottleneck effect. Despite moderate resolution, 
the inertial ranges converge, which is due to locality of 
hydrodynamic cascade in spacial scales, making it pos- 
sible to consider higher and lower resolution simulations 
on a common ground, neglecting the influence of large 
scales, where energy is provided by driving. 

Fig. [2] presents a resolution study for simulations Rl- 
3 determining the spectral slope and Kolmogorov con- 
stant for Alfvenic turbulence. If the spectrum —3/2 
was universal, the outer scale point, corresponding to 
k = 3 which is marked by a cross will go down from 
Rl to R3 by a factor of around 1.26, instead it stays 
at about the same level, indicating that deviations from 
—5/3 slope are small (note that the outer-scale point 
moves horizontally in Fig. Q] as well). The fiat part of 
Fig. [2] in R3 simulation between k = 54 {k-q ks 0.037) 
and k = 91 (kr/ « 0.063) with central frequency k — 70 
was fit to obtain Kolmogorov constant. The value ob- 
tained in this fit was Ck = 3.2 ± 0.2 where the error was 
mostly due to fluctuation of spectrum in time. 
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Fig. 3. — Dynamic alignment seems to saturate towards small 
scales before the dissipation scale. This indicates that the scale- 
dependent alignment is a transient effect that is present in simu- 
lations due to MHD turbulence being less local. The asymptotic 
regime of MHD turbulence is showing constant alignment, which 
will not modify the spectral slope —5/3 of strong MHD turbulence 
(GS95 slope). 

4. DYNAMIC ALIGNMENT 

It was suggested that the spectral slope of MHD turbu- 
lence is modified by so-called "dynamic alignment" that 
increases indefinitely towards small scales. Although the 
tentative correspondence with theoretical scaling from 
Boldyrev model has been obtained with only one par- 
ticular measure of alignment, this was interpreted by 
some studies as a confirmation of the aforementioned 
model. In this paper we refer to our earlier studies 
IBeresnvak h Lazarianl (2006, 2009b) that measured sev- 
eral types of alignment and their dependence on scale. 
In these studies there were no conclusive evidence that 
all alignment measures follow the same scaling. In this 
paper we confirm this finding with higher-resolution sim- 
ulations, in addition we found evidence that all align- 
ment measures saturate, i.e. approach an asymptotic 
constant value on small scales. Fig. [3] shows the 
alignment measures in R3, where AA, AA2, DA and 
PI are different alignment measures: AA = (|<5w^" x 
<5w-|/|Jw+||5w-|), AA2 = (\5v+ x ^|/|<Jv+||^|>, 
PI = (\Sw+ x^w-|)/(|<5w+||5w-|), DA = (\5v+ X 
<5b x |)/(|<5vtll£b % |), fo r more details and motivation, see 
IBeresnvak fc Lazarianl ()2006l . [2009b). Having an iner- 
tial range of around two orders of magnitude in scale, if 
align ment was proportional to A ' as in IBoldvrevI ( 2005. 
I2006T ). we would expect alignment increase by a factor of 
3.2, while in reality the polarization intermittency PI in- 
creases only by a factor of 1.3, and dynamic alignment 
DA by a factor of 1.8. This is consistent with the range 
of A between 3 and 10. 

We are not aware of any convincing physical argumen- 
tation e xplaining why al ignment should be a power-law 
of scale. IBoldvrevI (|2006l) provides an explanation argu- 
ing that alignment will tend to increase indefinitely, but 
will be bounded by field wandering, i.e. the alignment on 
each scale will be created independently of other scales 
(hence the term "dynamic alignment") and will be pro- 
portional to the relative perturbation amplitude 5B/B. 
But this directly violates precise symmetry of Eq. ® i.e. 
Eq. [3] which states that nothing should depend on SB/B 
as long as other quantities are scaled properly. Phys- 



4 



ically, this means that held wandering can not destroy 
alignment or imbalance. Indeed a perfectly aligned state, 
e.g., with (5w~ = is a precise solution of Eq. [2] and it 
is n ot destroyed by its own held w andering. Addition- 
ally, iBcrcsnva k fc Lazarianl ()2009bl ) measured alignment 
in simulations of strong MHD turbulence with different 
values of SB^/Bq and found very little or no dependence 
on this parameter. Fig. [3] also compares alignment mea- 
sure with a first-order structure function of t he pertur- 
bation amplitude, i.e. (|<5wa|). According to iBoldyrevi 
(2006) they should scale the same way, but this is not 
observed. 

To summarize, our numerical data are consistent with 
alignment measures becoming constant in the inertial 
range and inconsistent with the hypothesis that they de- 
pend as A 1 / 4 on scale. This finding is important, be- 
cause if alignment is constant on scale in the asymptotic 
regime, there is no reason to expect that the power-law 
scaling of turbulence will deviate from its —5/3 value 
for strong Goldreich-Sridhar turbulence. This result is 
further supported by the results of the previous section 
where a steeper asymptotic spectra has been observed. 

5. THE AMOUNT OF SLOW MODE AND THE TOTAL 
KOLMOGOROV CONSTANT FOR MHD TURBULENCE 

Full incompressible MHD turbulence have a cascade of 
slow mode, which was not included in our reduced MHD 
simulations Rl-3. Although in nature slow mode is of- 
ten damped, it is normally present in full MH D incom- 
pressi ble simulations, e.g. the ones presented in lBiskampI 
(2003). The passive cascade of slow mode will have the 
same energy spectral slope as a Alfvenic mode, and, as- 
suming that the ratio of slow to Alfvenic energies is C s , 
the total Kolmogorov constant for MHD turbulence will 
be expressed as 

C K = C KA {l + C s ) 1 '\ (5) 

The ratio C s is supposedly depend on how the MHD 
turbulence is driven. However, historically, previous 
studies simulated MHD turbulence with zero mean held, 
either d ecaying or driven with sta tistically isotropic forc- 
ing, e.g. iMuller fc Grappinl (|2005ft . In this idealized case 
Kolmogorov constant has been measured , although with 
fairly limited resolution (jBiskampl f2003). In this paper 
we will use a less straightforward approach, by measuring 
C s from a simulation with zero mean held and substitut- 
ing it into Eq. [5] This approach is motivated by our 
finding that MHD turbulence is less local and therefore 
it is much harder to achieve an asymptotic universal cas- 
cade if one uses zero-mean field simulation. Indeed, one 
has to observe the transition to the strong local mean 
field case, which will require at least a couple of order 
of magnitude in scale and subsequently a transition to 
universal cascade, which as we observed in previous sec- 
tion, takes about two orders of magnitude in scale, as 
long as the power-law scaling and Kolmogorov constants 
are concerned. It is, therefore, impossible to directly 
measure the properties of universal cascade in zero mean 
field simulations of currently available simulations. The 
"natural" value of C s is unity, because the incompress- 
ible MHD equations have four degrees of freedom, out of 
which Alfvenic mode uses two and slow mode also uses 
two. Having the same amount of degrees of freedom and 



the isotropic driving that does not prefer any direction we 
would expect that the energy will be distributed equally 
between the modes. We measured how energy is parti- 
tioned on small scales of simulation Ml by making a local 
Fourier transform of smaller cubes and decomposing into 
modes with respect to the local mean held. The actual 
partition of energy shows C s being around 1.3. Although 
statistical errors in this measurement are small, it is hard 
to claim a particular value of C s based on a numerical 
simulation with a finite resolution. Conservatively, we 
will assume that C s is between 1, which is equipartition, 
and 1.3, which is observed in our simulation Ml. The 
total Kolmogorov constant will be estimated as 4.1 ±0.3. 

6. SCALE LOCALITY AND KOLMOGOROV CONSTANT 

The energy flux through scales in both MHD and hy- 
drodynamic turbulence can be expressed as a certain 
third order signed structure function divided by scale and 
has to be scale-local due to an upper analytical bound 
on contributions fro m different k wavebands (see, e.g. 
lAluie fc Eivnkl 12010ft . This upper bound, however, is 
well applicable to similar third order unsigned structure 
function. This unsigned third order structure function 
is related by self-similarity hypothesis to second order 
structure function, which is a measure of energy. There- 
fore, we would expect that the ratio of unsigned third 
order structure function to the signed one will scale ap- 

3/2 

proximately as CJ . This seriously limits the bound on 
scale locality from practical standpoint as long as Ck 
beco mes large, i.e., the energy tran sfer becomes less effi- 
cient (|Beresnvak fc Lazarianll2010ft . Indeed, if the dchnc 
"scale locality" as a ratio of largest to smallest wavevec- 
tors k\jki which still significantly contribute to energy 
flux through some central wavevector ko, then this ratio 
will have an upper bound that scale asymptotically with 
Kolmogorov constant as C^ 4 . In other words, inefficient 
energy transfer can still be very local, but it is also pos- 
sible that it is less local than efficient energy transfer. 
A nonlocal or diffuse energy transfer must be inefficient 
and must have a high value of Kolmogorov constant. 

A comparative stu dy of energy spectra in MHD and 
hydro turbulence in iBeresnvak fc Lazarianl (|2009bft re- 
vealed that the bottleneck effect is less pronounced or 
altogether absent in MHD simulations, while in hydro 
it is always present, both in simulations with normal 
(n = 2) and hyperviscosity (n > 2). This was inter- 
preted as an indication that MHD cascade is less local. 
Now, our measurement of Kolmogorov constant revealed 
that MHD energy transfer is less efficient, therefore MHD 
cascade may be less local than hydro cascade. In view of 
all numerical evidence available today, MHD cascade is 
most likely less local than hydro cascade. 

7. DISCUSSION 

Previous measurements of the slope usually relied on 
the highest-resolution simulation and fitted the slope in 
the fixed fc-range close to driving scale typically between 
k = 5 and k = 20. In this paper we argue that such 
a fit is unphysical and instead one should fit a fixed krj 
range. In the former case the result would be a shallower 
spectral slope due to proximity to the outer scale and 
driving. In the latter case the effect of the driving will 
diminish with increasing resolution and one will observe 
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shallower spectra at small resolutions that will become 
steeper with increasing resolution. 

Earlier measurements of Kolmogorov constant in MHD 
turbulence re ported lower va lues than this study, e.g. 
C K = 2.2 in IBiskampI (|2003l ). We believe this is due 
to insufficient resolution in those simulations, which pre- 
vented the observation of the asymptotic regime. In par- 
ticular, in the case of s t atisti cally isotropic simulations 
like the ones in IBiskampI (|2003f) a transition to small scale 
subAlfvenic regime precede the transition to asymptotic 
regime. These two transitions require numerical resolu- 
tion that is even higher than the highest resolution pre- 
sented in this paper and for now seems computationally 
impossible. Our own statistically isotropic simulation 
Ml shows Kolmogorov constant of 3.5, which is still only 
a lower limit, consistent with 4.1 derived in this paper. 
For Ml and similar lower-resolution simulations the esti- 



mate of Ck continues to grow with increasing resolution, 
which supports argumentation above. 

In this paper we treated so called balanced case, 
where the rms amplitudes of the w ± were statistically 
the same. A number of attempts to generalize the 
GS95 model has been made recen t ly (iLithwick et all 
20071: iBeresnvak fe Lazarianl 120081 : iChandranl 120081: 
Perez fc Boldvrevl |2009D~ Some of these models can 
be rejected by numerics, in particular t he model based 
on alignment (jPerez fc Boldvrevl 120091 ) is grossly in- 
consistent with the dissipation rates measured in im- 
balanced nu merical simulations (jBeresnyak fc Lazarianl 
I2009al 1201(1 . 
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